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Using atomistic simulations we investigate the morphological properties of graphene deposited 
on top of a nanostructured substrate. Sinusoidally corrugated surfaces, steps, elongated trenches, 
one dimensional and cubic barriers, spherical bubbles, Gaussian bump and Gaussian depression are 
considered as support structures for graphene. The graphene-substrate interaction is governed by 
van der Waals forces and the profile of the graphene layer is determined by minimizing the energy 
using molecular dynamics simulations. Based on the obtained optimum configurations, we found 
that: (i) for graphene placed over sinusoidally corrugated substrates with corrugation wave lengths 
longer than 2nm, the graphene sheet follows the substrate pattern while for supported graphene 
it is always suspended across the peaks of the substrate, (ii) the conformation of graphene to the 
substrate topography is enhanced when increasing the energy parameter in the van der Waals model, 
(iii) the adhesion of graphene into the trenches depends on the width of the trench and on graphene's 
orientation, i.e. in contrast to a small width (3 nm) nanoribbon with armchair edges, the one with 
zig-zag edges follows the substrate profile, (iv) atomic scale graphene follows a Gaussian bump 
substrate but not the substrate with a Gaussian depression, and (v) the adhesion energy due to van 
der Waals interaction varies in the range [0.1-0.4] J/m^. 



INTRODUCTION 



Geometrically structured substrates affect various 
properties of graphene [l|, Q , and can prevent the crum- 
pling of graphene which is typical for free standing 
graphene without a support [3]. Before graphene's dis- 
covery in 2004, the study of 2D membranes over cor- 
rugated substrates was an important branch of soft- 
condensed matter physics with e.g. applications in bio- 
logical systems 0, [HI . Recently, particular attention was 
focused on the various properties of graphene on top of 
a substrate. Substrates can induce corrugations, modify 
the electric conductance and deform graphene [6, 7]. The 
electrostatic interaction of graphene on a substrate can 
be understood as due to the van der Waals (vdW) inter- 
action of graphene with the metallic gate below the sub- 
strate, the electrostatic forces between graphene and the 
polarized substrate, the water between the substrate and 
graphene, and impurities between graphene and the sub- 
strate [8]. The vdW interaction includes attractive and 
repulsive terms, which are widely used and extensively 
investigated in soft matter [9]. The usual dispersion in- 
teraction for the attractive part of the vdW interaction 
(sum of contributions proportional to D~^) must be mod- 
ified for the two 7r-conjugated systems at distance e.g. 
graphite [10]. For the vdW interaction between carbon 
nanostructures and Si-C/Si02 substrates the Lennard- 
Jones potential (LJ) is widely used and produced both 
qualitative and quantitative acceptable results pj]-[T6|. 
Ab-initio calculations obtained vdW energy curves be- 
tween carbon nanostructures (those curves are similar 
curves to the LJ function), e.g the vdW interaction be- 
tween hexagonal boron nitride sheets [17], methane ad- 
sorption on graphene and gas molecule adsorption in car- 
bon nanotubes were found to be LJ like functions when 
using vdW corrected density functional theory [l8|, [TqI . 



Recently, experimental measurements on the adhesion 
energy of pressurized mono-layer/multi-layer graphene 
on top of a Si02 substrate showed that the adhesion en- 
ergy is ultra strong (x ~ 0.3 — 0.45 J/m^) which is many 
times larger than the one reported for typical microme- 
chanical structures and is of the order of solid-liquid ad- 
hesion energies [20]. This adhesion energy is one order 
of magnitude larger than the upper limit found for water 
modified adhesion between graphene and the substrate. 
Kusminskiy et al used x=2 meVA~^ for the pinning 
of a tethered membrane (as a model for graphene within 
continuum elasticity theory) and found the possible mor- 
phology of graphene over a Gaussian bump and a Gaus- 
sian depression [21]. Their model includes both bending 
and stretching energies together with a constant pinning 
energy. 

Here, as distinct from previous works, we investigate 
graphene on top of several nanostructured substrates 
with different geometrical deformations. We carried out 
molecular dynamics simulations at T=300K to minimize 
the energy and find the optimum profile of the deposited 
graphene membrane. Sinusoidal substrates with differ- 
ent wave lengths, elongated trenches, barriers, bubbles, 
Gaussian bump and Gaussian depressions are considered 
as geometrical distinct examples of nanostructured sub- 
strates. We find that in case of a sinusoidal substrate 
with short wave length and small energy parameter in 
the vdW model (i.e. e), graphene does not follow the 
substrate. For graphene on top of the trench, it is found 
that zig-zag graphene falls into the well but arm-chair 
graphene is suspended across the trench. The stress dis- 
tribution shows that the atoms within the deformed parts 
are highly stressed. For the boundary conditions of the 
examined graphene flakes we considered both free and 
supported (i.e. fixed) in-plane boundaries. We found 
a significant difference in the obtained graphene profile 
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when on top of a Gaussian bump or at a Gaussian de- 
pression, i.e. the graphene sheet over a depression with 
1 nm variance and 1 nm height does not fah into the de- 
pression while it follows a Gaussian bump with the same 
size. For a Gaussian bump/depression with larger vari- 
ance, graphene follows both substrates. The square bar- 
rier (a cube with 1 nm side) influences graphene such that 
an unexpected pyramidal shape is found which surrounds 
the barrier. We studied the vdW energy stored between 
graphene and the nanoscale Gaussian bump by employ- 
ing a continuum model for both systems and calculated 
the variation of the vdW energy per area as a function 
of the energy parameter of the model. We also com- 
pared our molecular dynamics results for the Gaussian 
bump/depression to those predicted by the continuum 
model and found agreement only in case of a large Gaus- 
sian bump with weak interaction , i.e. small e. 

This paper is organized as follows. In Sec. II the de- 
tails of the atomistic model are presented. In Sec. Ill we 
present the continuum model for the vdW interaction of 
graphene and various substrates. In Sec. IV we present 
results for various nano-structured substrates and differ- 
ent boundary conditions. The results are summarized in 
Sec. V. 



II. ATOMISTIC MODEL 

Classical atomistic molecular dynamics simulation 
(MD) is employed to simulate large flakes of graphene 
(GE). The second generation of Brenner's bond-order po- 
tential is employed which is able to describe covalent sp^ 
bond breaking and the formation of associated changes in 
the atomic hybridization within a classical potential [22| . 
The Brenner potential terms were taken as 



(1) 
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where Ep is the average binding energy, and and 
are the repulsive and attractive terms, respectively, 
where Vij is the distance between the atoms i and j (all 
relevant parameters in this study are listed in Table I). 
Bij is called the bond order factor which includes all 
many-body effects. Bij depends on the local environ- 
ment of the bond, i.e. the bond and torsional angles, the 
bond lengths and the atomic coordination in the vicinity 
of the bond. This feature allows the Brenner potential 
to predict correctly the configurations, the hybridization 
and the energies for many different hydrocarbon struc- 
tures. 

The carbon-carbon bond length, ao, is 1.42 A. In our 
model, the origin of the Cartesian coordinate system 
is set at (0,0,0). Here, the two primitive vectors for the 
GE sublattices, ai = ^/^a^i and a2 = >/3/2aoi + 3/2aoj 
are the two basic vectors of GE lattice. 

In order to model the van der Waals (vdW) interac- 
tion between GE and different substrates, we employed 



the Lennard- Jones (LJ) potential. The LJ potential de- 
scribes both the repulsive and attractive parts of the vdW 
energy between two atoms which are non-bonded. The 
LJ potential is a widely used potential in various simu- 
lations [11-16J. For two interacting uncharged particles, 
we have 



u{r) = elaia/ry^ - f3{a/r) 



(2) 



where r is the distance between the two particles, e 
and a are the 'energy parameter' and the 'length pa- 
rameter\ respectively. Setting the integer numbers to 
a = P = m = 4 Eq. (2) gives the 12-6 LJ potential 
function and for a = 2, /3 = m = 3 Eq. (2) gives the 
9-6 LJ potential function [23]. The minimum of u{r) is 
Tmin = cr(^)i/(^^-^) which yields the 2^/^ and a for 
12-6 LJ and 9-6 LJ potential, respectively. Therefore, 
the equilibrium distance is shorter for the 9-6 LJ poten- 
tial while the minimum of u{r) {u\r^^^ = — e) is the same 
for both cases. Notice that for the 12-6 LJ potential both 
attractive and repulsive terms have the same weights in 
ix(r), i.e. a = p. We will use mostly the 12-6 LJ poten- 
tial (in some exceptional cases we use the 9-6 LJ potential 
will be mentioned expliciltly). 

To model the interaction between two different types of 
atoms such as the carbon atom (C) and a substrate atom 
(S), we adjust the LJ parameters using the equations 
= v'ece and ar = {crc + cr)/2. For carbon we use 
the parameters ac =3.369 A, and ec =2.63 meV. For the 
substrate atoms we vary a in the range (2.5A, 3.5A ) and 
e in the range (10.0 meV, 140.0 meV), where the lower 
limits are typically for insulators, e.g. Si02 [Hi and the 
upper limits are typical for metallic substrates, e.g. Na, 
K, etc [2^, |25|. Notice that the main difference between 
the two set of parameters is the energy parameter (e) 
which is varied over more than one order of magnitude. 

The total vdW energy stored between GE with N 
atoms and a substrate with M atoms, can be written 



as 



N M 



^vdW 



(3) 



where r^j 

■th 



\ri — Tjl and refers to the position of the 
carbon atom of GE and rj refers to the j^^ atom of the 
substrate. Often in MD simulations, one approximates 
the above sums by including only the nearest neighbor 
atoms in order to reduce the number of interactions. Such 
an approach is accurate in the case of short-range poten- 
tials. Regarding the cutoff distance appropriate for the 
LJ potential Tc = 3cr, only those substrate atoms inside 
a sphere having radius Tc with origin at the position of 
the i^^ atom of GE, interact most strongly with the i^^ 
atom, while outside this sphere, the interaction strength 
decreases very fast. Therefore, in practice for each 'i\ the 
sum over M can be truncated and limited to the atoms 
inside a sphere with radius Tc. This is done by employ- 
ing a neighbor list in our MD simulation. In our study. 
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TABLE I: List of all relevant parameters used in the paper. 





The length and width of the simulated graphene membrane 


e,cr 


The energy and length parameters in the van der Waals (vdW) potential for the substrate atoms, Eq. (|2]) 


a, /3, m 


Integer numbers defining the power law of repulsive and attractive parts of the vdW potential in Eq. (|2]) 




The stress and the force on atom i 


/, EydW 


Total free energy (Eq. (|6])) and vdW energy (Eqs. (|3)).(p^) stored in graphene when on top of a substrate 




Mass, volume, velocity and position of atom i 


N,M 


Total number of atoms respectively, of the simulated graphene and the substrate 


r, hv 


Surface tension, bending rigidity of the graphene membrane 


X,T 


The vdW energy per planar area 'A' and temperature 


hG{x,y),hs{x,y) 


Height of the graphene membrane and the substrate at (x,y) 


Ep,B^j,VR,VA 


Total bond energy of graphen, bond order, repulsive and attractive terms in the Brenner potential, Eq. ([1]) 




Hamaker constant, the density of the simulated substrate and graphene membrane, respectively 


X,0{x) 


The wave length and the step function 




The radius of Gaussian bump and the determinant of the metric tensor 


ho 


The amplitude of the sinusoidal waves or height (depth) of Gaussian bump/bubble/barrier (depression or trench) 


hi, d 


A shift or vertical distance between graphene and the substrate and the width of the trenches/barriers 


, Train 


The Lennard- Jones potential energy and its minimum distance 


ao, ai, a2 


Carbon-carbon bond length and two basic vectors of the graphene lattice 


I 


Lattice parameter of the substrate lattice 


rc,ro 


Cutoff distance in vdW interaction, upper limit of the integrals in Eq. ()10p 




Distance between a lattice site of graphene ('i') and a lattice site of the substrate ('j') 


Ari2 


Distance between a surface element of the graphene membrane and a surface element of the substrate 



the number of GE atoms is N =14400 (which is equiv- 
alent to a graphene sheet with dimension Ix =19.17nm 
and ly =19.67 nm) and the number of substrate atoms 
is M =6000 (only in the case of the Gaussian bump we 
performed a simulation with A/'=72000 and M=35000). 

The adhesion energy can be obtained using ab-initio 
calculations and can be estimated using classical models 
(e.g. LJ potential). The present day patterns of de- 
formation of large scale GE on top of substrates is be- 
yond reach of traditional ab-initio methods. Our clas- 
sical vdW model as based on the LJ potential are able 
to simulate the realistic sizes and gives a vdW energy 
(the main term in the binding energy) between two non- 
bonded systems. Note that several ab-initio calculations 
have demonstrated that the vdW interaction between 
nanoscale objects can be well approximated by a LJ po- 
tential [l3-[i9j. 

The gradient of the total potential energy, i.e. Etotai — 
Ep + E^ind^ force experienced by the i^^ carbon 

atom, = —ViEtotai • In common molecular dynamics 
simulations Newton's second law should be solved nu- 
merically in order to determine the path of motion of 
the atoms. In this study, the equations of motion were 
integrated using a velocity- Verlet algorithm with a time 
step of 0.5 fs and the temperature was held constant at 
T=10K by a Nose-Hoover thermostat. 

The atomic stress experienced by each i^^ atom can be 
expressed as [26|, [l^ 

vU = ^(^\rnvlvi + ^rr^Ft}j, (4) 

where the inner summation is over all the carbon atoms 
which are neighbors of the i^^ atom which occupies a vol- 



ume Q = 47rao/3. The quantities m and v'^ denote the 
mass and velocity of the i^^ atom and the scaler r^j is 
the u component of the distance between atoms 'i' and 
'j'. Fj^j is the force on the i^^ atom due to the j^^ atom in 
the /i direction. We have used this expression to calculate 
the stress on each atom. In order to be able to visual- 
ize the stress distribution on the GE atoms, we colored 
the atoms according to the value of the dimensionless 
stresses, i.e green (red) is related to the minimum (max- 
imum) possible stress. 

III. CONTINUUM APPROACH 

Evaluation of the vdW contribution of the stored en- 
ergy in the deformed GE with average density Sg, due 
to the interaction with a substrate with average density 
E^, is obtained after the integration of the vdW poten- 
tial over both GE and the substrate surfaces. Here, we 
present for comparative purposes a continuum model for 
the stored vdW energy between the GE membrane and 
the substrate. Such an approach can be used to calculatie 
the vdW energy stored between two objects 0, 

In the absence of external pressure, the total free en- 
ergy of a membrane consists of three terms, i.e. bending, 
stretching and vdW terms which are given by 

f= JJ dxdy[T{VhG{x,y)f ^ n{V^hG{x,y)f]^ E^^^, 

(5) 

where r and n are the surface tension and the bend- 
ing rigidity of the membrane, respectively. The two first 
terms in Eq. ([5]) are relevant to the bending and stretch- 
ing energies of GE [21] and the third term, E^^^^^ is the 
total vdW contribution of the interaction between the 
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substrate and the membrane. E^^^^ includes two repul- 
sive and attractive terms. The stored adhesion energy 
per area is determined by [4] 



XT = ( 



fm 



A 



r), 



(6) 



where /^^^ is the minimum of the total free energy when 
the membrane takes its optimum configuration. 'A' is 
the projected area onto the x — y plane, i.e. A = J dxdy. 
Equation (|6|) was used by Swain et al [J, [sj to estimate 
the adhesion energy of typical soft membranes over dif- 
ferent substrates. In Ref. [2l| the adhesion part in the 
free energy was taken as a coupling constant. Assum- 
ing a planar local relative height coordinate function for 
the vdW interaction energy between the substrate and 
the soft membrane, i.e. the Deryagin approximation, the 
vdW energy is approximated by 



e: 



c 

vdW 



= 11 dxdy[Vo 



(7) 



where ip is di function of the height increment Sh = 
hcix^y) — hs{x^y) and Vq is a constant. Substituting 
Eq. ([7j) in Eq. (j5j) and minimizing the total free energy 
with respect to haix^y) results in the following differen- 
tial equation |5j 



{kV"^ - rV^ + v)hG{x, y) = vhs{x, y), 



(8) 



where v is proportional to the Hamaker constant (cx 
ecr^SsSq). Equation ^ can be solved in Fourier 
space Ulil 



/iG(k) 



/is(k) 



(9) 



where k is the wave vector, = r/v and = hz/v. 
Equation ([9]) was used to find the optimum configuration 
of a soft membrane on top of corrugated substrates |5| . 

Here we assume that the vdW energy is not local- 
ized and use the LJ potential for the interaction between 
graphene and the substrate. This gives us the vdW con- 
tribution to the adhesion energy, i.e. x- We assume that 
both substrate and membrane are homogenous and con- 
tinuous materials. The obvious difference between the 
atomistic model and the continuum model is the absence 
of the chirality effect. The LJ potential energy between 
the GE membrane and the continuum substrate is given 
by 



7C 

^vdW 
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C(ri)C(r2)ii(Ari2) dsids2, 



(10) 



where Ari2 = |ri - r2| and C(ri) = SG(ri)^G(ri) and 
C(r2) = S5(r2)^5(r2). Here, Eg = 2/|ai x a2| and 
are the mean surface density of carbon atoms in the GE 
and the substrate lattice, respectively. The substrate 
density is = £~'^ which is equivalent to a (100) sur- 
face of a crystal with lattice parameter equal to £. The 



scaler ^(r^) = y 1 + h{xi^yi)Y is defined by the ap- 
propriate transformation of a surface element in curvi- 
linear coordinates into a two-dimensional planar surface 
{x — y plane in cartesian coordinate), i.e. the determi- 
nant of the metric tensor . In Eq. (p!Q|) the position vector 
= {xi^yi^ h{xi^yi)) refers to the position of the surface 
element dsi of the GE membrane {i = 1) or substrate 
{i = 2). As already mentioned, because of the short 
range nature of the LJ potential, we assume that the 
main contribution of the vdW energy is due to the in- 
teraction with the outer surface of the substrate. This 
model gives a good insight into the vdW adhesion en- 
ergy between GE and various substrates. For the contin- 
uum model we will give only the results for GEs on top 
of a Gaussian bump/depression. Notice that it is ana- 
lytically impossible to minimize Eq. (|5|) by substituting 
E^^^^ (Eq. (do])). Therefore, we used the optimum con- 
figuration for the GE membrane as obtained from our 
MD simulations. 



IV. RESULTS AND DISCUSSION 

In this study we investigate several different geome- 
tries for the substrate which can be realized experimen- 
tally. The substrate atoms are assumed to be rigid during 
the simulations which is a reasonable approximation due 
to the different atomic- vibrations time scale in graphene 
and the substrate. In order to model the substrate, 
a (100) surface having lattice parameter ^=3 A is used 
which is a typical lattice parameter. The density of sites 
in the substrate is = £~'^. Since the interaction be- 
tween the substrate and graphene is weak and of short 
range (i.e. van der Waals interaction) the main contribu- 
tion to the mutual stored energy and to the force between 
graphene and the substrate comes from the upper layer 
of the substrate. We can show this explicitly by calcu- 
lating the energy and force as function of the number of 
atomic layers in the substrate for one of our samples (We 
took the case of graphene over a Gaussian bump see Fig. 
12(a)). The total energy can be written as 



L N,M 



(11) 



n=l i,j 



where Un is the contribution of the n^^ layer and L is the 
number of considered substrate layers. For instance for 
graphene on top of a Gaussian bump the attractive and 
repulsive parts are proportional to (p^ + (z + n^)^)~^ and 
(p^ + (2; + n£)^)~^, respectively (here p {z) is the planar 
(vertical) distance between an atom in graphene with one 
in the top layer of the substrate) which decrease fast with 
n. Recalling the discussion below Eq. (3) we found that 
Eq. (fTTj) for the considered system around the central 
points (r <30 A) the top layer (n = 1) contributes almost 
99% of the total energy, the second layer contributes 1% 
and the contribution of the other layers are neglect ible. 



FIG. 1: (Color online) Side views of the optimum configura- 
tion of graphene on top of sinusoidal substrates with different 
wave lengths. The filled white circles are the substrate atoms. 
The parameters in (a,b,c) are a — 3.5 A, e = lO.^meV and 
in (d,e,f) are cr=3.4A and e = 100.0 meF. The wave lengths 
are A =2nm (a,d), A =3nm (b,e) and A =4nm (c,f). For the 
graphene sheet, the colors represent the stress distribution, 
the highest stress is denoted by red and the lowest by green. 



This motivated us to restrict our study to the top sub- 
strate layer which helps us considerably to minimize the 
CPU time. Note that for larger and smaller cft the 
contribution of the second layer will increase. 

The height of the graphene and the substrate atoms at 
each point {x^y) are denoted by hcix^y) and hsix^y)^ 
respectively. The calculations are done for two different 
boundary conditions: i) free boundary condition, and ii) 
supported boundary conditions which prevents in-plane 
movements for two longitudinal ends of GE. The two 
atom rows at the longitudinal ends were taken in the 
zig-zag direction (in most cases) and they are allowed to 
move in the z-direction. 



A. Sinusoidal substrates: Free boundary condition 

A sinusoidal deformation of the substrate along the 
x-direction is given by 



hs{x,y) = hosm{kx), 



(12) 



where k = 27r/A and the amplitude is h^. We used dif- 
ferent wave lengths, i.e. A =2, 3 and 4nm and two sets 
of cr and e, i.e. (3.5 A, 10.0 meV) and (3.4 A, 100.0 meV) 
with fixed ho = 0.5 nm. At the start of the simulation, 
we put a flat graphene sheet on top of this substrate at 



hcix.y) = hs{x,y)- 



i. We choose the x-direction to 



be the arm-chair direction. 

After 0.5 ns of the MD simulation, GE found its op- 
timum configuration which is deformed and corrugated. 
Figure [1] shows six snap shots of free GE (upper cor- 
rugated sheets in each panel) over three different sub- 
strates (filled circles in each snap shot). In Figs. [Ha,b,c) 



FIG. 2: (Color online) Two side views of Fig. H^c). For 
the graphene sheet, the colors indicate the stress distribution. 
The substrate atoms are indicated by filled white circles below 
graphene. The sinusoidal substrate is shown only by the first 
front row of atoms. 



the vdW parameters were set to cr=3.5A, e=10.0meV 
and in Figs. [Hd,e,f) the vdW parameters were set to 
cr=3.4A, 6=100.0 meV. The wave length in Figs. ma,d) 
is 2 nm. Notice that for A=2 nm GE does not follow the 
substrate, i.e. hcix^y) ^ hi ^ hs{x^y)^ where hi is a 
vertical shift of the order of graphite's layer distance, 
i.e. 3.4 A. By increasing the wave length or e, GE fol- 
lows more closely the substrate profile, i.e. hcix^y) = 
hi + hs{x^y). Increasing e yields stronger adhesion and 
deforms GE (Fig. [2] shows zoomed versions of Fig. ttl^c)). 
Therefore, according to our MD simulations, the short- 
est wave length which makes GE's profile similar to the 
substrate's profile is larger than A =2nm. In all figures 
the colors indicate the stress distribution. Notice that at 
the boundaries we always have red colors (i.e. maximum 
stress) because of the presence of dangling bonds. In the 
other parts, the distribution of stress is uniform (green 
colors). Notice that in each particular system we scaled 
the colors by the highest stress. 

Here we compare our results to those predicted by con- 
tinuum elasticity theory for a membrane on top of sinu- 
soidal surfaces. The possible solution of Eq. (j8j) [4] for a 
membrane on top of a deformation given by Eq. ([12]) is 



hcix.y) 



ho sm{kx) 



(13) 



For longer wave lengths, or small and r this solution 
gives hcix^y) ^ hs{x^y) which is in agreement with 
our MD results for long wave lengths and large e. For 
graphene membrane k, ^ 1.1 e V and r ~ eV/A^ (taken 
of the order of the Lame coefficients) and Eq. (p!3]) is 
valid if V n^r which are related to both large e, 
and E^. The stiffer membrane with larger n and r can 
not curve easily and stronger adhesion due to a larger v 
(oc e, i.e. stronger adhesion) is required. By using a = 
[crc + crs)/2 = (3.369 + 3.4)/2 A, e = V2.63 x 10 meV, 
= 0.225 A"^ and = 0.026 A-^ yields the Hamaker 
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FIG. 3: (Color online) The optimum configuration of 
graphene over a substrate with wave length A = 4nm. 
Graphene is supported by the longitudinal ends while it can 
move along the ^-direction. Here e = lOmeV and a = 
3AA for the 12-6 LJ potential (see also Fig. [D^c) which shows 
graphene with free boundaries over the same substrate). 
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constant ~ Air'^ea^TiG^s = 1-78 eV and consequently 
V 0.1 and thus k'^^r < 1 and k^^^ <C 1 yield A > 4.7 A 
and A ^ 5. 6 A, respectively, which are in agreement with 
our MD finding, i.e. A > 20A. 



B. Sinusoidal substrates: Supported boundary 
condition 



Here, we impose the supported boundary condition on 
the longitudinal ends (which is mostly taken to be the 
zig-zag direction). In this case the atoms at x = 
are not allowed to move in-plane while they are al- 
lowed to move in the z-direction. In this section we use 
the 12-6 LJ potential with parameters e =10meV and 
a =3.4 A. We repeated the above simulations by applying 
the supported boundary condition along the zig-zag di- 
rection for graphene on top of a substrate with A =4 nm. 
Fig. [3] shows the optimum configuration of the deposited 
graphene over this substrate after minimization. Notice 
that the obtained deformation is different from those for 
free graphene over the same substrate. Fig. [D^c). The 
reason is the supported boundary (i.e. fixed in-plane) at 
the two longitudinal ends. Graphene does not follow the 
substrate, but the lateral edges (along the x-direction at 
±/^/2) feel a much large stress as compared to Fig.lTJc). 
Graphene is suspended across the periodic peaks with a 
small curvature between them. The larger e enhances the 
latter effect. Therefore, the vdW energy is not dominat- 
ing the bending energy of GE. 





FIG. 4: (Color online) Arm-chair (a,b) and zig-zag (c,d) 
graphene on top of two steps of height 1 nm shown along two 
different angle of view. The highest stressed atoms are shown 
by red and the lowest by green. The substrate atoms are 
shown by filled circles below graphene. 



Figure m shows two snap shots of the optimum configu- 
rations (an arm-chair GE with two different directions of 
view in (a,b) and a zig-zag GE with two different direc- 
tion of view in (c,d)) of free GEs on top of steps. All GEs 
follow the steps except around x ^ where GE is bent in 
a continuous fashion and does not follow the substrate. 
There are no considerable differences between the opti- 
mum configurations of the free arm-chair (Figs. Ill^a,b)) 
and free zig-zag (Figs. |4l^c,d)) GEs on top of the step. 
Indeed, the wall atoms at x = shift (due to adhesion) 
in both the left and the right parts of the GE towards 
x = 0. 



C. Step: Free boundary condition 

The second class of interesting substrate configurations 
are steps which were recently studied in an experiment 
to measure the electronic and morphology of deposited 
graphene [29| 

hs{x,y) = hoO{x), (14) 

where 6{x) is the Heaviside step function with step height 
of ho =1 nm and with density of sites. Both GE with 
arm-chair and zig-zag direction are put on top of the 
steps. 



D. Step: Supported boundary condition 

In Fig. [5] the optimum configuration of GE along the 
arm-chair direction with supported boundary condition 
is shown which is over a sharp step defined by Eq. (fT4|) . 
Notice that there is a significant difference between the 
deformation obtained here and the one depicted in Fig|4j 
The curvature around the step {x ^ 0) are different and 
all atoms of GE feel more or less equal stress. The wall 
at X = adheres both the left and the right part of the 
GE but the supported ends prevent fully adhesion of GE 
to the wall, especially for the right hand side of GE. 
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FIG. 5: (Color online) The optimum configuration of arm- 
chair graphene over a step located at x = with sup- 
ported longitudinal edges (see also Figs. ma,b) which show 
graphene with free boundaries over the same substrate) . Here 
e =10.0 meV and a =3.4 A for the 12-6 LJ potential. 



E. 



Trench: Free boundary condition 



The other important substrate that we studied here is 
an elongated trench 



hs{x,y) = hoO{x^ - d?), 



(15) 



with two walls located at x = id = ±1.5 nm with step 
height of 1 nm and with T^s density of sites. Figure [6] 
shows two snap shots of GE on top of such trenches (an 
arm-chair GE with two different angle of view in (a,b,c) 
and a zig-zag GE with two different angle of view in 
(d,e,f)). After MD minimization zig-zag GE follows the 
trench except around x ~ ±lnm. In this region, zig-zag 
GE is bent and does not follow the substrate. There is 
a significant difference between the optimum configura- 
tions of arm-chair (Figs.[6l^a,b)) and zig-zag (Figs.[6l^d,e)) 
GEs. An arm-chair GE does not follow the trench as well 
as a zig-zag GE. We attribute this effect to the larger 
number of atoms of zig-zag GE (as compared to arm- 
chair GE) in the well region (\x\ < d), which yields a 
strong attractive force on the GE atoms due to the sub- 
strate atoms within the well's wall. Recently, Lu et al [23] 
studied wider trenches with 2d =28.6 nm using a 9-6 LJ 
potential in order to find the vdW adhesion of GE mem- 
branes. In our study we also used the 9-6 LJ potential 
and found different deformations as compared to the 12-6 
LJ potential, see Figs. [6fc,f). This is due to the different 
strength of both attractive and repulsive parts in the two 
models. 



F. Trench: Supported boundary conditions 

In Fig. [71 we show the optimum configuration of arm- 
chair (a) and zig-zag (b) graphene with supported bound- 
ary condition on top of the trench defined by Eq. (fT5]) . 
There is a significant difference between the deformation 
obtained here and those depicted in Fig. [6l The curva- 
ture for |x| < d is very different and GE atoms around the 
well feel a lower stress as compared to the one shown in 
Fig.O Here both arm-chair and zig-zag GE do not follow 
the substrate and were suspended over the wells which is 
a consequence of the supported boundaries. Therefore, 




Arm-Chair 12-6 LJ 
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FIG. 6: (Color online) Arm-chair (a,b,c) and zig-zag (d,e,f) 
graphene on top of trenches 1 nm deep and 3 nm wide along 
the ^/-direction. Panels (a),(b) (and also (d),(e)) are the same 
with different angle of view. For a graphene sheet, the colors 
indicate the stress distribution. The substrate atoms are indi- 
cated by filled circles below the GE. Zig-zag graphene follows 
the trench in contrast to arm-chair graphene. In (a,b,d,e) and 
(c,f) we used the 12-6 LJ and the 9-6 LJ potential, respec- 
tively. 



by controlling the boundary condition one can clearly 
control the GE deformation over the substrate. 



G. Barriers: Free boundary condition 

A barrier in the middle of the substrate is the inverse 
situation of the previous ones. An elongated barrier in 
the ^/-direction is parameterized as 

hs{x,y) = h^e{x^ -d^), (16) 

with two walls at x = id = ±1.5 nm with step height 
of 1 nm and with T^s density of sites. Figs. [Ua,b) shows 
two snap shots of arm-chair GE (with two different angle 
of view) on top of the elongated barrier. As we see the 
stressed regions are located around x = ^ d. GE does 
not follow the rectangular shape of the barrier in this 
part. 

Another interesting case is the one of a substrate that 
consists of a cubic barrier in the middle (see the inset (i) 
in Fig. He)) 



hs{x,y) = h^e{x^-d^)e{y^-d^), 



(17) 
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FIG. 7: (Color online) The optimum configuration of arm- 
chair (a) and zig-zag (b) graphene over two trenches located 
at |x| < 1.5 nm where both ends were supported in the x — y 
plane (see also Figs. [6ja,d) which shows two graphene mem- 
branes with free boundaries over the same substrate). Here 
e=10meV and a =3.4 A for 12-6 LJ potential. 



with four walls at y = id = ±1 nm with step height of 
1 nm and with density of sites. Figs.[8fc,d) show that 
the optimum configuration is pyramidal shaped (inset (ii) 
shows this schematically). This particular deformation is 
due to the four corners of the cube which strongly repel 
the GE.. The highest stresses are distributed around the 
steps (red colors in the + ?/^| < 3d). The C-C bond 
lengths in GE are distributed non-uniformly (Fig. [Sfe)) 
but still around the barrier we have a larger stretch in 
the bond lengths (up to 0.147 nm which are shown by 
red colors). 



H. Barriers: Supported boundary conditions 

Fig. [9] shows the optimum configuration of arm-chair 
GE in the case of supported boundary condition over 
two different barriers (defined by Eqs. ([16]) , (fT7|) ). As 
we see the stress distribution and the deformations are 
completely different to those shown for free graphene (see 
Fig [5]). 



I. Spherical bubble: Free boundary condition 

The next important type of deformation for the sub- 
strate that has been realized experimentally [30|, is a 
bubble (see Fig. [TOT b)) which we model by 




/itiitm\\\^ 




FIG. 8: (Color online) Elongated barrier (a,b) and a cubic 
barrier (c,d) deformation in the middle of a substrate covered 
with graphene. The colors indicate the stress distribution in 
graphene. The inset (i) in (c) shows a schematic model for the 
cubic barrier and the inset (ii) is a schematic for the optimum 
configuration of graphene, i.e. pyramidal shape. Panel (d) is 
another view of (c). The C-C bond lengths distribution for 
(c) (or (d)) are shown in (e), where the red colors are related 
to the bonds around 1.47A and green colors are related to the 
bonds around 0.140A. Here e = lOmeV and a = 3AA for 
12-6 LJ potential. 
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FIG. 9: (Color online) The optimum configuration of arm- 
chair graphene over an elongated barrier of size |x| < 1.5 nm 
where zig-zag edges were supported in the x — y plane (see 
also Fig. [8] which shows graphene with free boundaries over 
the same substrate). Here e = lOmeF and a = 3.4A for 12-6 
LJ potential. 



hs{x, y) = ^JB? - + hx, 



(18) 



where R is the radius of the bubble and = ^y^. In 
order to create an uniform discrete atomistic structure for 
the bubble, we set hi = —R/2. Figures [TOl (ax) show the 
obtained optimum configuration from MD simulation for 
the GE on top of a bubble with R=2 nm. The optimum 
configuration is a Gaussian. In order to produce uniform 



bubbles, we increased the density of lattice sites in the 
bubble location where S = 1/4 Increasing the num- 
ber of lattice sites on the bubble results in a much larger 
stressed GE which influences regions far from the center 
(see Fig.IlOl^a)). 
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FIG. 10: (Color online) (a) The optimum configuration of the 
graphene sheet over a bubble with hs{x,y) = y^R^ - p^-R/2 
deformation (R = 2nm). The corresponding substrate is 
shown in (b). In (c) we show both the GE and the sub- 
strate. For a graphene sheet, the colors indicate the stress 
distribution. The highest stress is shown by red color and the 
lowest by green. The substrate atoms are indicated by filled 
white circles below graphene. 




FIG. 11: (Color online) The optimum configuration of arm- 
chair graphene with two longitudinal ends supported in the 
X — y plane on top of a bubble (see also Fig. [10] which 
shows graphene with free boundaries over the same sub- 
strate). The inset shows the elongation of the deformation 
of graphene along the x-direction, i.e. arm-chair direction. 
Here e — lOmeV and a =3.4A for 12-6 LJ potential. 



J. Spherical bubble on the substrate: Supported 
boundary condition 

The optimum configuration for the supported 
graphene over a bubble substrate is shown in Fig. [TTJ 
Due to the supported ends GE elongates longitudinally 
along the supported direction, see inset in Fig. [TTJ 

K. Gaussian bump/depression: Free boundary 
conditions 

There have been already a few studies that evaluated 
different properties of a GE membrane in the presence 
of a Gaussian deformation, but those studies did not ad- 
dress the following issues: (i) the creation of the Gaus- 



sian deformation in GE using an atomistic scale deformed 
substrate; (ii) what is the effects of the vdW energy 
strength on both the deformation and the adhesion en- 
ergy at the atomistic scale; (iii) what are the important 
differences between the deformation due to a bump and 
due to a depression on the atomistic scale, and (iv) what 
is the effect of the boundary condition on GE. 

The Gaussian bump (protrusion) /depression [U, [32| 
are parameterized as (Fig. [T2fb)) 

hs{x, y) = ±ho exp(-pV27'), (19) 

where -\-ho{—ho) is the height (depth) of the Gaussian 
bump (depression) and = + y'^ and 7 is the vari- 
ance of the Gaussian. Kusminskiy et al studied recently 
the pinning of GE over a Gaussian bump in order to find 
the corresponding attachment /deattachment of GE [21]. 
Our model is more realistic with relevant length scales 
for both height and variance of the bumps/depressions. 
Figure [EJ^a) shows a snap shot of the optimum con- 
figuration of a GE on top of the Gaussian bump (see 
Fig. [T2r b)) with /iq = 7 = 1 nm located at the center 
where e =10.0 meV and a =3.4A. The inset of Fig.[T2]^a) 
shows the far view of GE and the inset of Fig.[T2fb) shows 
the side view of the Gaussian bump. The red colors refer 
to the highest stresses which are mostly located around 
the bump region, r < 27. For this size of the bump, GE 
follows the Gaussian bump, i.e. hs{x^ y) ^ hi-\- hcix, y), 
where hi is a vertical shift which is about graphite's layer 

spacmg 3.4A. 

Figure \T2\c) shows the optimum configuration of GE 
(found from MD) on top of a depression with ho=-l nm, 
and 7 = 1 nm. As we see the deformation of the GE over 
the Gaussian bump is different from the one over the de- 
pression (while both have the same variance, heights and 
potential parameters, i.e. e =10.0 meV and a =3.4A). 
This is clear from the curves shown in Fig. [13] which 
were taken along x = and ^ = (corresponding to the 
deformations shown in Figs.[T2ja,c)). The optimum con- 
figuration of GE on top of the depression (Figs. [THTcd)) 
is not a Gaussian, i.e. hsix^y) ^ hi ^ hcix^y)^ because 
of the stronger repulsive force inside the depression due 
to the interaction between GE and the substrate. 

Figure [15] shows the variation of both E^^y^ (MD) 
(Eq. (O) and E^^y^ (CM) (Eq. ([TOj)) versus the radius, 
r, where r is the upper limit of the integrals in Eq. ([TO]) . 
In Fig. [T5lf a). we set 7 = /io=lnm, e = lOmeV and 
a = 3. 4 A which are close to the one for the Si02 sub- 
strate [11]. For the substrate with a bump the energies of 
the atomistic model (MD) are close to the one obtained 
from the continuum model (CM). For the depression, 
our MD results are different from the continuum model 
results, which is a consequence of the different profiles 
in GE and the substrate (see Fig. [T3]) . Notice that the 
used profile in CM for graphene on top of a depression is 
a Gaussian profile (compatible to the substrate profile) 
while the found optimized profile in the MD simulation 
as seen from Figs. 13(c,d) is not a Gaussian. Therefore, 
there is a significant deviation from the CM and MD re- 
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FIG. 12: (Color online) The optimum configurations obtained 
from molecular dynamics simulations for a graphene sheet 
places on top of a Gaussian bump (a) Gaussian depression 
(c) with hs{x,y) = 10 exp(-pV200)(A). Here we see that 
graphene on top of the Gaussian bump follows the substrate 
which is not the case for a depression. For a graphene sheet, 
the colors indicate the stress distribution. Here e = lOmeF 
and cr = 3.4A. Panel (b) is the bumped substrate, i.e. 
hs{x,y) = 10exp(V/200)(A). 




FIG. 13: (Color online) The heights of the substrate (red) 
and graphene (blue) along x = and ^ = on top of a Gaus- 
sian bump/Gaussian depression in x-direction (a) /(c) and y- 
direction (b)/(d) (they are taken from the central portion of 
those deformations shown in Fig. [12)). Circular points indi- 
cate the C-atoms of graphene and the solid curves are the 
substrate. Notice that graphene follows the bump which is 
not the case for a depression. The colors indicate the stress 
distribution in the graphene sheet,. Here e =10.0 meV and 
cr =3.4A. 




FIG. 14: (Color online) Bond length distribution in graphene 
on top of a Gaussian bump (a) and a Gaussian depression 
(b). The longer bond lengths are shown by red colors and the 
shorter bond lengths are shown by green colors. 



suits for depression as shown in Fig. 15(a). In the inset of 
Fig.fTsTa). which is related to graphene on top of a bump, 
we set the energy parameter to e = 140 meV. When the 
energy parameter is large, the results of MD and CM de- 
viate, which is related to the strong attraction between 
the substrate and GE. 

Figure [TST b) shows the vdW energy for a Gaussian 
bump/depression with larger variance, i.e. 7 =3nm and 
/io=lnm, e =10meV and cr =3.4A. There is good agree- 
ment between the results of MD simulations and those 
found from the continuum model (CM). We conclude 
that for large variance the continuum model provides 
a vdW contribution to the adhesion energy which are 
comparable to the MD atomistic results. But for small 
bump/depression the CM model is not applicable and 
the lattice structure of GE should be taken into account. 

In Fig. [16] the variations of x = ^vdw/^ versus e for 
various a (=2.5 A, 3.5 A) are shown. Here, graphene is 
deposited on top of a Gaussian bump with 7=1 nm, and 
/io=lnm where r = ro =2nm and the area is calculated 
using A = Trrg = 47rnm^. In Fig. [ToT a) and Fig. [ToT b) 
we used the 12-6 LJ and the 9-6 LJ potential parameters, 
respectively. The 9-6 LJ potential gives results that have 
larger |x| for a particular e (Notice that in this paper the 
9-6 LJ potential is used only for comparative purposes). 

The energy per area, i.e. Xi is in the range of the adhe- 
sion energy found for a graphene membrane positioned on 
top of Si02 substrate [20], i.e. 0.31-0.45 J/m^ However, 
note that our results give only the vdW contribution of 
the adhesion energy which results from the standard 
dispersion interaction. 



L. Gaussian bump/depression: Supported 
boundary conditions 

Since the optimum configuration of supported 
graphene over a Gaussian bump is similar to the one 
for a spherical bubble, we will not report them here. For 
supported graphene over a Gaussian depression the opti- 
mum configuration is not Gaussian, similarly as for free 
graphene over a depression (we do not show the optimum 
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r(A) r(A) 

FIG. 15: (Color online) (a) Variation of the van der Waals energy versus the radius (in log-scale) measured from the origin 
(0,0). Results are presented for both molecular dynamics simulation (MD, Eq. (|3])) and continuum model (CM, Eq. ()10p ). 
Here the energy parameter in the LJ potential was set to 10 meV (in the inset we took 140 meV) and a = SAA. The substrate 
is a Gaussian bump with 10exp(— /> V200)(A) deformation, (b) Variation of the vdW energy versus the radius (r) for wider 
Gaussian bump/depression with 10 exp(— p^/900)(A) deformation with e =10meV and a — 3.4A. 




^(meV) E(meV) 



FIG. 16: (Color online) Variation of van der Waals energy per area versus e from continuum model (Eq. ()10p ) for graphene 
on top of a Gaussian bump with 10 exp(— p^/200) deformation, where a increases from top to bottom with steps of 0.25 A. In 
(a),(b) we used 12-6 and 9-6 LJ potential parameters, respectively. In both panels r = ro = 2nm which is the upper limit of 
the integrals in Eq. ([TO)) (see the text). 



configuration here). V. SUMMARY 

We carried out several molecular dynamics simulations 
and studied systematically the optimum configuration 
of large scale graphene deposited on top of several 
differently shaped substrates. The stress distribution in 
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graphene shows that highly stressed atoms are located 
around the deforraed regions of the substrate. 

For short wave length (< 2nm) graphene is suspended 
across the neighbor peaks of the sinusoidal substrate and 
thus graphene does not follow the substrate. A graphene 
sheet on top of a cubic barrier shows an unexpected 
pyramidal shape. It is found that for large Gaussian 
bump/depression the van der Waals contribution in the 
adhesion energy are in agreement with the prediction 
of the continuum model. The van der Waals adhesion 



energy per area for a nanoscale Gaussian bump is found 
to be around 0.1-0.35 J/m^ depending on the energy 
parameter of the model. 
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